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Abstract 

Two dimensional simulations of plane and ax- 
isymmetric jets are presented. These simula- 
tions were made by solving full Navier-Stokes 
eqautions using a high order finite difference 
scheme. Simulation results are in good agree- 
ment with the linear theory predictions of the 
growth of instability waves. 

1. Introduction 

The recent efforts to develop accurate nu- 
merical schemes for transition and turbulent 
flows are motivated, among other factors, by 
the need for accurate prediction of flow noise. 
The success of developing high speed civil trans- 
port plane (HSCT) is contingent upon our un- 
derstanding and suppression of the jet exhaust 
noise. The radiated sound can be directly 
obtained by solving the full (time-dependent) 
compressible Navier-Stokes equations. How- 
ever, this requires computational storage that 
is beyond currently available machines. This 
difficulty can be overcome by limiting the solu- 
tion domain to the near field where the jet is 
nonlinear and then use acoustic analogy [e.g., 
Lighthill 1 ] to relate the far-field noise to the 
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near-field sources. The later requires obtaining 
the time-dependent flow field. 

The other difficulty in aeroacoustics compu- 
tations is that at high Reynolds numbers the 
turbulent flow has a large range of scales. Di- 
rect numerical simulations (DNS) cannot ob- 
tain all the scales of motion at high Reynolds 
number of technological interest. However, 
it is believed that the large scale structure 
is more efficient than the small-scale struc- 
ture in radiating noise [see, e.g., Bishop et. 
al 2 , Crington 3 , Goldstein 4 , Hussain 5,6 , Kibens 7 , 
Liu 8,9 , Mankbadi and Liu 10 , Mankbadi 11,12 , 
Mollo-Christensen 13 , and Zaman 14,15 ]. Thus, 
one can model the small scales and calculate 
the acoustically active scales. The large scale 
structure in the noise-producing initial region of 
the jet can be viewed as a wavelike nature, the 
net radiated sound is the net cancellation after 
integration over space. As such, aeroacoustics 
computations are highly sensitive to errors in - 
computing the sound sources. It is therefore es- 
sential to use a high-order numerical scheme to 
predict the flow field. 

The present paper presents the first step in a 
ongoing effort to predict jet noise. The empha- 
sis here is in accurate prediction of the unsteady 
flow field. We solve the full time-dependent 
Navier-Stokes equations by a high order finite 
difference method. Time accurate spatial simu- 
lations of both plane and axisymmetric jet are 
presented. Jet Mach numbers of 1.5 and 2.1 
are considered. Reynolds number in the sim- 
ulations was about a million. Our numerical 
model is based on the 2-4 scheme by Gottlieb 
& Turkel 16 * 17 . Bayliss et al. 18 applied the 2-4 
scheme in boundary layer computations. This 
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scheme was also used by Ragab and Sheen 19 to 
study the nonlinear development of supersonic 
instability waves in a mixing layer. 

In this study, we present two dimensional di- 
rect simulation results for both plane and ax- 
isymmetric jets. These results are compared 
with linear theory predictions. These computa- 
tions were made for near nozzle exit region and 
velocity in span wise /azimuthal direction was as- 
sumed to be zero. Computational domain of the 
present study is shown in Figure 1. 

2. Governing Equations 

The flow field of the technologically impor- 
tant high-Reynolds-number compressible jet is 
governed by the compressible Navier- Stokes 
equations, which can be written as 


LQ = 0 (1) 

Where L is a three-dimensional operator and Q 
is the solution vector. Here we present the two 
dimensional form of the operator L pertaining 
to our formulations in the polar coordinates. 
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3. Numerical Scheme 
3.1 Discretization 

For problems in computational acoustics one 
needs to use high order (at least fourth or- 
der) schemes. Hence, in this study, Gottlieb 
and TurkeFs extension 16 * 17 of MacCormack’s 
scheme is used. The extension is fourth-order 
accurate in space and second order accurate in 
time and is known as the 2-4 scheme. For three 
dimensional computations, operator L in equa- 
tion (1) is split into three one-dimensional op- 
erators and the 2-4 scheme is applied to each of 
the three operators as 

Predictor step 

Q t = Q7 + ~m F ,\i-K)~ 

(f? +a - (8) 

Corrector step 


<?i n+1 = \iQi + Qi + W - *T-i) 

-(f?*- 1 - F?-2)} + A<S,] (9) 

This scheme becomes fourth-order accurate 
in the spatial derivatives when alternated with 
symmetric variants. We define L\ as an one 
dimensional 2-4 scheme operator with forward 
difference in the predictor and backward differ- 
ence in corrector. Its symmetric variant Li uses 
backward difference in predictor and forward 
difference in the corrector. 

For two dimensional computations, one di- 
mensional sweeps are arranged with alternate 
symmetric variants as 

<? n+1 = L lx L lr Q n (10.1) 

Q n+J = L 2r L 2x Q n+1 (10.2) 

Our scheme has a truncation error of the 
form 0(At((Ax) 4 + (At) 2 + (A<)(Ax) 2 )). For 
A t — 0((Ax) 3 ) as Ax — * 0, this scheme be- 
comes fourth order accurate. 


F and G are the fluxes in x and r directions 
respectively, and S is the source term that arises 
in the cylindrical polar coordinates, and r,j are 
shear stresses. 


3.2 Boundary Conditions 
The 2-4 scheme uses one sided differences of 
the fluxes. For a computational domain extend- 
ing from i=l to m, fluxes are needed at m-hl, 
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m+ 2 , - 2 , >1 in addition to the interior points. 
The fluxes at points outside the computational 
domain are estimated using cubic extrapolation 
from the interior, i.e., 

Fm +1—4 F m — 6 F m _i + 4F m _3 - F m _ 3 ( 11 . 1 ) 

-Fm+2 — 4T m +i — 6T m — F m _2 (11*2) 

At the centerline, a new set of equations 
are derived from the original equations using 
L’Hospit&ls rule to circumvent numerical prob- 
lems associated with geometric singularity in 
the formulation. 

Physical boundary conditions for the com- 
putations are derived using linearized charac- 
teristics. For supersonic flows, all character- 
istics travel in the flow direction. At inflow 
all variables are given. For outflow, the vari- 
ables are calculated applying the 2-4 scheme at 
the boundary. Extrapolations of variables out- 
side the domain are done using equation ( 11 ). 
We stress that extrapolation of fluxes to artifi- 
cial points is identical to using one sided differ- 
ences. The extrapolation is used for program- 
ming convenience only. For subsonic flows, one 
characteristic variable propagates against the 
flow direction while the rest follow the flow di- 
rection. For inflow, three characteristic vari- 
ables are specified and the other one is extrap- 
olated from interior. We formulated our out- 
flow boundary conditions following the results 
of Bayliss and Turkel 20 . In particular, we solve 
equation ( 12 ) to predict the flow variable at the 
boundaries. 


P t - pcu t = 71 

( 12 . 1 ) 

P t + pcu t = *12 

( 12 . 2 ) 

Pt - c 2 pt = 73 

(12.3) 

v t - 74 

(12.4) 


The derivatives of P y u, v y p are then con- 
verted to derivatives of the conservative vari- 
ables. The right hand side of the above equa- 
tions are calculated from the solution obtained 
by applying the 2-4 scheme at the boundaries. 
If the flow at a point at the outflow boundary 
is subsonic, yi is set to zero. If the flow at the 
boundary is supersonic, the value of is kept 


unchanged. Equation (12) is then solved to get 
corrected temporal derivatives at the bound- 
ary. Thus, for supersonic flow ( 12 ) is equiva- 
lent to using the PDE at the outflow with one 
sided differences. For subsonic flow a nonre- 
flecting boundary condition is used. In a fu- 
ture paper several nonreflecting boundary con- 
ditions will be compared. These all are in- 
cluded in generalizations of ( 12 . 1 ). We ob- 
tained the present form of ( 12 . 1 ) by simplifying 
the Bayliss-Turkel formulation 30 and neglecting 
spatial derivatives. Flow at the top boundary 
is always subsonic and a similar characteristic 
boundary condition is applied. 

4. Results 

In this study, our primary goal is to simu- 
late aerodynamic noise associated with a super- 
sonic jet flowing into a subsonic free stream. In 
the short term, we are interested in comput- 
ing growth rates of the disturbances imposed at 
the inflow. These computations are done in two 
stages, (a) first a mean (steady) state of the 
field with steady inflow condition is obtained 
and then (b) transient behavior of the flow with 
periodic excitations at the inflow is calculated. 

For the mean state calculations, we start with 
an initial field which is homogeneous in the axial 
direction. The initial axial velocity is specified 
as 

®(y) = '[(1+ «oo) - (1 - Uoo) 

tanh(a(y — yi )] (13) 

and corresponding temperature is specified by 
the Busemann-Crocco integral of energy equa- 
tion as 

T(y) = + (To - TooK* - u^/i 1 - u TO ) 

+.5(7 - 1)M 3 (1 - w)(u - tx*,) (14) 

The parameter M a” in the equation (13) deter- 
mines the sharpness of the velocity profile. The 
axial velocity is normalized by its initial value 
at the center of the jet. At location yi at the 
inflow plane, axial velocity is the average of the 
jet center and free stream value. Throughout 
the whole simulation, variables at the boundary, 
including those at the inflow are updated using 
characteristics boundary condition described in 
§ 3.2. We assume Prandlt number to be 1 , and 
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the kinematic viscosity is calculated using Su- 
therland law. The initial static pressure is as- 
sumed uniform across the field and density is 
calculated from the equation of the state. Ini- 
tial radial/ transverse velocity in the field is set 
to be zero. 

Once the steady state of the field is reached, 
a time varying disturbance is applied at the in- 
flow. In our present study, only the axial veloc- 
ity is perturbed. It is done by prescribing the 
axial velocity at the inflow as 

u(y,t) = u(y)(l -f € $in(u?t)) (14) 

Other variables at the inflow are kept at their 
steady state values. Depending on the value 
of the forcing frequency (a?), input disturbance 
undergoes either amplification or decay. For an 
amplifying disturbance, the initially growth is 
linear. However, as the amplitude of the distur- 
bance grow, the process becomes nonlinear and 
disturbances with other frequencies appear in 
the field. Smaller values of e delays the appear- 
ance of nonlinear effects downstream. Corre- 
sponding to a forcing frequency (u>) of the input 
disturbances, there are some disturbance eigen- 
functions ($ w ) of the flow variables. Distur- 
bances with shapes of the eigenfunction exhibit 
dominant growth. The shape of the disturbance 
as given in equation (14), is different from those 
of the eigenfunctions. With our choice of input 
disturbance, we left the simulation to generate 
the correct mode shape of the instability eigen- 
functions. Study of instability mode growth 
with input disturbances as eigenfunction shapes 
in all the variables are in progress. In the follow- 
ing we present simulation results in this study, 

4.1 Plane Jet 

A case of a coflowing plane jet is considered 
first. The jet Mach number, based on axial 
velocity at the center of jet was 1.5, velocity 
ratio (uoo/uq) and temperature ratio (Too/T 0 ) 
were .74 and 2 respectively. This combination 
made Mach number at the free stream about 
the half of the jet Mach number. Parameter 
“a” [equation (13)] in the initial axial velocity 
profile was 4. The computational domain ex- 
tended 50 jet thickness in the downstream direc- 
tion and 2.5 jet thickness in the transverse direc- 
tion. We used 600 mesh points in axial and 60 
mesh points in transverse direction. The mesh 
was uniform in axial and stretched in transverse 
direction. Reynolds number based on the jet 


thickness and the inlet axial velocity was 1.27 
million. 

4.1.1 steady state simulation 

Parameters for our simulation were chosen 
to given very small spreading of the jet. This 
condition is ideal for comparison of the simula- 
tion results with weakly nonparallel linear the- 
ory. Profiles of axial velocities at three down- 
stream locations are shown in Figure 2. These 
velocities were computed after the steady state 
was reached. Differences in these profiles are 
insignificant. Contour plots of steady state vor- 
ticity and axial velocity for the whole domain is 
shown in Figure 3. The y axis in these figure is 
magnified to show the the whole domain. The 
flow remains virtually parallel and spreading of 
the jet or any tread towards the formation of 
the potential core is not noticeable. If the flow 
is not well resolved, numerical (truncation) er- 
rors inject disturbances in the flow and that in 
turn could give rise to different flow structures 
in the solution. Our choice of grid results in 
well resolved flow for the present study. Small 
linear growth of both momentum and vortic- 
ity thicknesses was observed in the steady state 
flow field. Vorticity thickness as defined by 

S v = (U 0 - Uoo)/(du/ dy) max 


varied between .5 and .505. Also, momentum 
thickness defined as 


f pjUp - U){U - gg) 

J (u 0 - u^y y 


varied between .086 and .087. Once the steady 
state solution is established, we concentrate on 
the unsteady state simulation. 

4.1.1 Time dependent simulation 

Now we shall discuss results of a time de- 
pendent simulation with periodic disturbance 
for w — tt/ 5 and e = .001. Contour plots of 
vorticity and axial velocity are shown in Fig- 
ure 4. Input disturbance grows spatially and 
causes vorticities to form. These structures be- 
comes prominent downstream. Near the out- 
flow, the input disturbance reached its maxi- 
mum strength and made the flow to oscillate. 
Care should be taken in boundary treatment, 
since numerical reflections at the boundary also 
contribute to such oscillations. 


Figure 5 shows mean axial velocities at three 
downstream locations. These were obtained by 
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computing the zeroth Fourier mode of the un- 
steady axial velocity. As expected, mean axial 
velocity remained at their steady state values. 
After long time, the flow field became spatially 
nonhomogeneous but time periodic about its 
steady state. Since the shape of the inflow per- 
turbation did not correspond to any eigenfunc- 
tion of instability wave, the flow underwent an 
adjustment region where the simulation picked 
up the correct mode shape of instability waves 
downstream. In Figure 6 the mode shape of 
the eigenfunction (<£) corresponding to the forc- 
ing frequency (u>) at two downstream locations 
are compared with the linear theory prediction. 
Mode shapes form the direct numerical simula- 
tions are obtained by taking Fourier transform 
in time. As expected, DNS prediction of these 
shapes improved as the flow moved away from 
the inlet region. There is, however, a small 
phase shift between the DNS prediction and the 
linear theory predictions. In Figure 7, we ex- 
amine the growth of the instability mode ($ w ) 
corresponding to the forcing frequency, whose 
shape is given in Figure 6. For the DNS, growth 
rate is calculated along y=t/i line. At this verti- 
cal location, the shear in the mean axial velocity 
profile is the highest. There is a good agreement 
between the simulation and the linear theory. 
DNS predictions differ from the linear theory 
predictions near outflow. This may be due to 
excitation of instability waves at other frequen- 
cies or non-linear effects in the flow. Growth 
rate predicted by both the linear theory, 

and DNS at X =50 is .054. As the distur- 
bances grow spatially downstream, non-linear 
effects excited other instability modes. Two 
such modes for the simulation are shown in Fig- 
ure 8. For our simulation, any such mode was 
significant only near the outflow boundary. 

Simulations were also performed for Jet Mach 
number 2.1 and higher velocity ratios. Com- 
parisons of one such simulation results with the 
linear theory are shown in Figures 9 fe 10. Ve- 
locity and temperature ratios were .2 and 1 re- 
spectively. Parameter w a” in the initial axial 
velocity profile was 6. Computational domain 
extended 35 jet thickness downstream and 2.5 
jet thickness in the transverse direction. Mesh 
size in the axial and the transverse directions 
were 400 and 150, and stretching was used in 
both directions. For the unsteady calculations, 
we used w — x/8 and e = 10~ 5 . 


4.1 Axisymmetric Jet 

Results of the axisymmetric jet simulations 
are similar to those in the plane jet case. Here 
we present results of an axisymmetric jet simu- 
lation with Mach number 1.5 in the jet, veloc- 
ity ratio (uqo/uo) = .75 and temperature ratio 
/To) = 2. “a” in the axial velocity pro- 

file equation was 4. The computation domain 
extended 100 radii downstream and 5 radii in 
radial direction. 400 grid points were used in 
axial and 100 grid points were used in radial di- 
rection. The grid spacing was uniform in axial 
and stretched in radial direction. The Reynolds 
number of this flow was the same as in plane jet 
case. 

4-1.1 steady state simulation 

Contour plot of the vorticity for the steady 
state simulation is shown in Figure 11. The 
Y (radial) axis in these figures is magnified to 
show the whole computational domain. As in 
the case of plane jet, in steady state simulation, 
the flow remained virtually parallel. 

4.1.1 Time dependent simulation 

A time dependent simulation with w = x/4 
and e — .005 was made to study the growth of 
disturbance mode. Contour plots of vorticity 
and axial velocity of the time dependent field 
are shown in Figure 12. These plots show spa- 
tial growth of oscillatory flow structures caused 
by instability modes. Growth of disturbance 
corresponding to the forcing mode is compared 
with the linear theory predictions in Figure 13. 
Agreement between the linear theory and DNS 
predictions are similar to those in plane jet case. 
Small discrepancy between linear theory and 
DNS near inlet is likely due to the fact that 
the flow underwent an adjustment region, be- 
cause the input disturbance did not have correct 
shape of an instability eigenfunction. Discrep- 
ancy near outflow region are likely because of 
the similar reasons as in plane jet case. Growth 
rate (-a;) predicted by both the linear theory 
and DNS at X =50 is .035. Comparison of mode 
shape as predicted by the linear theory and DNS 
simulations at 50 radii downstream is given in 
Figure 14. Except close to the center of the jet, 
DNS prediction of the mode shape is in good 
agreement with the linear theory. Excitation of 
two modes due to nonlinear effects are shown in 
Figure 15. DNS results in Figures 14 & 15 are 
obtained exactly as their counterparts in plane 
jet simulations. 
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5. Conclusions 

In this study, we have presented a set of 
Naviei Stokes simulations for both two dimen- 
sional plane and axisymmetric jets. Parame- 
ters were so chosen to give almost parallel mean 
flow. This enabled us to compare our results 
with weakly nonparallel linear theory. Agree- 
ment with the linear theory are quite good. We 
imposed time periodic disturbances which did 
not correspond to the instability eigenfunction 
shape. Nevertheless, the simulation was able to 
generate the correct mode shape after an adjust- 
ment region (~ 10 diameters) and comparison 
with linear theory afterward is good. 

The boundary condition is very important 
for jet simulation. We formulated the outflow 
boundary condition by simplifying the Bayliss- 
Turkel 20 nonreflecting boundary condition. Ef- 
forts to improve the outflow boundary condition 
are in progress. 
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